library(tidyverse)
library(here)
library(phyloseq)
library(vegan)
theme_set(theme_bw())
filt_fun <- function (x, min_reads = 100, min_samples = 5) {
(sum(x) > min_reads) & (sum(x > 0) > min_samples)
}
ps <- readRDS(here('data','following_study','ps_rarefied.rds'))
sample_data(ps)[is.na(sample_data(ps))] <- 'NA'
# transform data into proportion
ps.prop <- ps %>% transform_sample_counts(function(x) x/sum(x))
# double check proportion transform
sample_sums(ps.prop) %>% summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
sam <- data.frame(sample_data(ps))
max.core <- parallel::detectCores()
alpha.diversity <- data.frame(sample_data(ps),estimate_richness(ps))
ggplot(alpha.diversity,aes(x = as.numeric(Household), y = Shannon, group = Household)) +
geom_point() + geom_line() + xlab('Household')
anova(lm(Shannon~Household, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.18604 1.5916 0.05425 .
## Residuals 49 5.7277 0.11689
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity,aes(x = Epileptic.or.Control, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
Paired t test
epileptic <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
control <- alpha.diversity %>% filter(Epileptic.or.Control == 'Control')
t.test(epileptic$Shannon, control$Shannon, paired = TRUE)
##
## Paired t-test
##
## data: epileptic$Shannon and control$Shannon
## t = -0.22089, df = 48, p-value = 0.8261
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.1556564 0.1248409
## sample estimates:
## mean difference
## -0.01540773
ANOVA with Household effect added
anova(lm(Shannon~Household + Epileptic.or.Control, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.186045 1.5607 0.06329 .
## Epileptic.or.Control 1 0.0058 0.005816 0.0488 0.82612
## Residuals 48 5.7219 0.119206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity) +
geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
facet_wrap(~Epileptic.or.Control) +
theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())
# here breed group with na is removed
anova(lm(Shannon~Household + Breed.Group..1., data = alpha.diversity %>% filter(Breed.Group..1. != 'NA')))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 45 7.4046 0.16455 1.5646 0.07819 .
## Breed.Group..1. 4 1.2267 0.30668 2.9160 0.03339 *
## Residuals 39 4.1016 0.10517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alpha.diversity.epi <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(alpha.diversity.epi, aes(x = Pheno.Y.N, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
t.test(Shannon~Pheno.Y.N, data = alpha.diversity.epi)
##
## Welch Two Sample t-test
##
## data: Shannon by Pheno.Y.N
## t = 0.98407, df = 33.649, p-value = 0.3321
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
## -0.1316909 0.3787797
## sample estimates:
## mean in group No mean in group Yes
## 2.856864 2.733320
ggplot(alpha.diversity, aes(x = Sex, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova(lm(Shannon~Household + Sex, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.9301 0.18604 1.6309 0.04673 *
## Sex 1 0.2522 0.25218 2.2107 0.14360
## Residuals 48 5.4755 0.11407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(Shannon~Sex, data = alpha.diversity)
##
## Welch Two Sample t-test
##
## data: Shannon by Sex
## t = -0.83555, df = 88.411, p-value = 0.4057
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -0.2223635 0.0907199
## sample estimates:
## mean in group F mean in group M
## 2.792318 2.858140
ggplot(alpha.diversity, aes(x = as.numeric(Age..months.)/12, y = Shannon)) +
geom_point() + xlab('Age')
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
lm(Shannon~as.numeric(Age..months.), data = alpha.diversity) %>% summary()
## Warning in eval(predvars, data, env): NAs introduced by coercion
##
## Call:
## lm(formula = Shannon ~ as.numeric(Age..months.), data = alpha.diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90685 -0.30999 0.01342 0.28351 0.75221
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.745594 0.091949 29.860 <2e-16 ***
## as.numeric(Age..months.) 0.000838 0.001114 0.752 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3838 on 94 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.005986, Adjusted R-squared: -0.004588
## F-statistic: 0.5661 on 1 and 94 DF, p-value: 0.4537
plot_ord <- function(data, factor, method, distance){
data.ord <- ordinate(data, method = method, distance = distance)
p <- plot_ordination(data, data.ord, color = factor)
p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
model <- as.formula(paste0('dist.matrix~', formula))
if (!is_null(strata)) {strata <- df[,strata]}
result <- adonis2(model,
data = df,
permutations=permutations,
strata = strata,
parallel = core)
return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
beta.disp <- betadisper(dist.matrix, group = df[,group])
result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
return(result)
}
permanova(ps.prop, 'Household', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.2596 9.999e-05 ***
## Residual 49 5.5192 0.31119
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_ord(ps.prop, 'Epileptic.or.Control','MDS','bray')
plot_ord(ps.prop, 'Epileptic.or.Control','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2235 -- in the phylogenetic tree in the data you provided.
plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','bray')
## Run 0 stress 0.2086595
## Run 1 stress 0.2041723
## ... New best solution
## ... Procrustes: rmse 0.05570532 max resid 0.4115003
## Run 2 stress 0.2049336
## Run 3 stress 0.2224345
## Run 4 stress 0.2036331
## ... New best solution
## ... Procrustes: rmse 0.03306726 max resid 0.305915
## Run 5 stress 0.204184
## Run 6 stress 0.2031278
## ... New best solution
## ... Procrustes: rmse 0.01942513 max resid 0.0845663
## Run 7 stress 0.2098507
## Run 8 stress 0.2099703
## Run 9 stress 0.2217262
## Run 10 stress 0.2258699
## Run 11 stress 0.2040622
## Run 12 stress 0.2229979
## Run 13 stress 0.2110726
## Run 14 stress 0.2042057
## Run 15 stress 0.2264617
## Run 16 stress 0.229214
## Run 17 stress 0.2194904
## Run 18 stress 0.2031184
## ... New best solution
## ... Procrustes: rmse 0.01466582 max resid 0.08599131
## Run 19 stress 0.2053133
## Run 20 stress 0.2155131
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 8: no. of iterations >= maxit
## 12: stress ratio > sratmax
plot_ord(ps.prop, 'Epileptic.or.Control','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV75 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1827303
## Run 1 stress 0.1948716
## Run 2 stress 0.1826111
## ... New best solution
## ... Procrustes: rmse 0.02702614 max resid 0.2334849
## Run 3 stress 0.1901819
## Run 4 stress 0.1827303
## ... Procrustes: rmse 0.027013 max resid 0.2336186
## Run 5 stress 0.1878658
## Run 6 stress 0.1882704
## Run 7 stress 0.1910674
## Run 8 stress 0.1821596
## ... New best solution
## ... Procrustes: rmse 0.03709371 max resid 0.2358942
## Run 9 stress 0.1813226
## ... New best solution
## ... Procrustes: rmse 0.04120045 max resid 0.2334379
## Run 10 stress 0.1904447
## Run 11 stress 0.1870262
## Run 12 stress 0.1879809
## Run 13 stress 0.1869531
## Run 14 stress 0.1809398
## ... New best solution
## ... Procrustes: rmse 0.04006268 max resid 0.2347613
## Run 15 stress 0.1896303
## Run 16 stress 0.1883852
## Run 17 stress 0.1896765
## Run 18 stress 0.1966149
## Run 19 stress 0.189808
## Run 20 stress 0.1831356
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.prop, 'Epileptic.or.Control', 'bray', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.1554 0.00876 0.8487 0.1039
## Residual 96 17.5803 0.99124
## Total 97 17.7358 1.00000
permanova(ps.prop, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2981 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.02437 0.01471 1.4336 0.05139 .
## Residual 96 1.63226 0.98529
## Total 97 1.65664 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.prop, 'Epileptic.or.Control', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00167 0.0016733 0.1117 10000 0.7396
## Residuals 96 1.43851 0.0149844
permdisp(ps.prop, 'Epileptic.or.Control', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1119 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0002514 0.00025138 2.0869 10000 0.1522
## Residuals 96 0.0115640 0.00012046
a table of # of breed. frequency table
(data.frame(sample_data(ps.prop))$Breed.Group..1.) %>% table()
## .
## Herder NA Pointer Spaniel Retriever Scent Hound
## 20 9 26 26 5
## Sight Hound Sled Dog Terrier
## 3 6 3
plot_ord(ps.prop, 'Breed.Group..1.','MDS','bray')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2583 -- in the phylogenetic tree in the data you provided.
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','NMDS','bray')
## Run 0 stress 0.2086595
## Run 1 stress 0.2055186
## ... New best solution
## ... Procrustes: rmse 0.05947923 max resid 0.4053172
## Run 2 stress 0.2217881
## Run 3 stress 0.2086136
## Run 4 stress 0.2223882
## Run 5 stress 0.2042025
## ... New best solution
## ... Procrustes: rmse 0.01751769 max resid 0.1099105
## Run 6 stress 0.2041345
## ... New best solution
## ... Procrustes: rmse 0.01516255 max resid 0.08075781
## Run 7 stress 0.2111371
## Run 8 stress 0.218612
## Run 9 stress 0.2030952
## ... New best solution
## ... Procrustes: rmse 0.03875405 max resid 0.3146495
## Run 10 stress 0.2121138
## Run 11 stress 0.2043533
## Run 12 stress 0.2041224
## Run 13 stress 0.2134545
## Run 14 stress 0.2041461
## Run 15 stress 0.2260212
## Run 16 stress 0.2087215
## Run 17 stress 0.2249867
## Run 18 stress 0.2053203
## Run 19 stress 0.2050026
## Run 20 stress 0.2046693
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 3: no. of iterations >= maxit
## 17: stress ratio > sratmax
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.prop, 'Breed.Group..1.','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV731 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1910295
## Run 1 stress 0.1947317
## Run 2 stress 0.1939297
## Run 3 stress 0.1892023
## ... New best solution
## ... Procrustes: rmse 0.04324271 max resid 0.2330391
## Run 4 stress 0.1915156
## Run 5 stress 0.1943251
## Run 6 stress 0.1884852
## ... New best solution
## ... Procrustes: rmse 0.03014957 max resid 0.2137321
## Run 7 stress 0.1919418
## Run 8 stress 0.1943292
## Run 9 stress 0.1877478
## ... New best solution
## ... Procrustes: rmse 0.06971062 max resid 0.2875466
## Run 10 stress 0.1928138
## Run 11 stress 0.1900625
## Run 12 stress 0.1932147
## Run 13 stress 0.1910529
## Run 14 stress 0.195628
## Run 15 stress 0.1929301
## Run 16 stress 0.1967131
## Run 17 stress 0.1891158
## Run 18 stress 0.1968316
## Run 19 stress 0.1999649
## Run 20 stress 0.1943933
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 18: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
# same breed group v.s. diff breed group in household
# check if data are paired within household
if (!identical(sample_data(ps)$Household[seq(1,98,2)],
sample_data(ps)$Household[seq(1,98,2)+1])) stop('data is not paired')
same.breed.group <- table(sample_data(ps)$Household, sample_data(ps)$Breed.Group..1.) %>%
apply(1, function(x) any(x == 2)) # if dog's grouping are same with in each house
same.breed.group;sum(same.breed.group)
## 1 10 11 12 13 14 15 16 17 18 19 2 20
## FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
## 21 22 23 24 25 26 27 28 29 3 30 31 32
## TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE
## 33 34 35 36 38 39 4 40 41 42 43 44 45
## TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 47 48 49 5 50 51 6 7 8 9
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [1] 39
Here we see 39 out of 49 household have dogs in same breed group
permanova(ps.prop, 'Household + Breed.Group..1.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.3397 9.999e-05 ***
## Breed.Group..1. 6 0.8417 0.04746 1.2896 0.07309 .
## Residual 43 4.6775 0.26374
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1666 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 0.82865 0.67272 2.2283 9.999e-05 ***
## Breed.Group..1. 6 0.07000 0.05683 1.5059 0.05389 .
## Residual 43 0.33314 0.27045
## Total 97 1.23179 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Here we
ps.control <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Control')
ps.epileptic <- ps.prop %>% subset_samples(Epileptic.or.Control == 'Epileptic')
perm.breed.bray.ctl <- permanova(ps.control, 'Breed.Group..1.', 'bray')
perm.breed.wunif.ctl <- permanova(ps.control, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2158 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.ctl;perm.breed.wunif.ctl
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.7534 0.194 1.4098 0.0295 *
## Residual 41 7.2851 0.806
## Total 48 9.0385 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.0033166 0.21306 1.5858 0.0386 *
## Residual 41 0.0122498 0.78694
## Total 48 0.0155664 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm.breed.bray.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'bray')
perm.breed.wunif.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2069 -- in the phylogenetic tree in the data you provided.
perm.breed.bray.epi;perm.breed.wunif.epi
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.3363 0.15645 1.0863 0.26
## Residual 41 7.2055 0.84355
## Total 48 8.5418 1.00000
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.21498 0.17043 1.2033 0.178
## Residual 41 1.04644 0.82957
## Total 48 1.26142 1.00000
\[X^2=-2 \sum_{i=1}^k \ln \left(p_i\right)\]
bray.p <- c(perm.breed.bray.ctl$`Pr(>F)`[1], perm.breed.bray.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(bray.p)), df = 4, lower.tail = FALSE)
## [1] 0.04501879
wunif.p <- c(perm.breed.wunif.ctl$`Pr(>F)`[1], perm.breed.wunif.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(wunif.p)), df = 4, lower.tail = FALSE)
## [1] 0.0410838
With using the weighted-unifrac distance, there’s at
least one of the null hypothesis is rejected.
Here we are only focusing on epileptic dogs.
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','bray')
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2547 -- in the phylogenetic tree in the data you provided.
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','bray')
## Run 0 stress 0.2020935
## Run 1 stress 0.1992304
## ... New best solution
## ... Procrustes: rmse 0.04145127 max resid 0.1509698
## Run 2 stress 0.2182018
## Run 3 stress 0.1994009
## ... Procrustes: rmse 0.06905381 max resid 0.3691135
## Run 4 stress 0.2034522
## Run 5 stress 0.2002274
## Run 6 stress 0.2117299
## Run 7 stress 0.2094952
## Run 8 stress 0.2029557
## Run 9 stress 0.2082424
## Run 10 stress 0.2087795
## Run 11 stress 0.2081331
## Run 12 stress 0.2041039
## Run 13 stress 0.2025866
## Run 14 stress 0.2106487
## Run 15 stress 0.2086674
## Run 16 stress 0.2145508
## Run 17 stress 0.2136685
## Run 18 stress 0.2061055
## Run 19 stress 0.2084818
## Run 20 stress 0.2100444
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 3: no. of iterations >= maxit
## 17: stress ratio > sratmax
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV3027 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1803735
## Run 1 stress 0.18135
## Run 2 stress 0.1940171
## Run 3 stress 0.1831955
## Run 4 stress 0.192244
## Run 5 stress 0.181868
## Run 6 stress 0.1814495
## Run 7 stress 0.1856995
## Run 8 stress 0.1827014
## Run 9 stress 0.1859298
## Run 10 stress 0.1803533
## ... New best solution
## ... Procrustes: rmse 0.097442 max resid 0.3798744
## Run 11 stress 0.1785945
## ... New best solution
## ... Procrustes: rmse 0.1067785 max resid 0.5529312
## Run 12 stress 0.1873117
## Run 13 stress 0.1856336
## Run 14 stress 0.186715
## Run 15 stress 0.1780467
## ... New best solution
## ... Procrustes: rmse 0.09321676 max resid 0.2595204
## Run 16 stress 0.1857093
## Run 17 stress 0.1849503
## Run 18 stress 0.1945162
## Run 19 stress 0.1873817
## Run 20 stress 0.1782206
## ... Procrustes: rmse 0.0993426 max resid 0.270433
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.epileptic, 'Pheno.Y.N','bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.2211 0.02588 1.2488 0.183
## Residual 47 8.3207 0.97412
## Total 48 8.5418 1.00000
permanova(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV794 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.01552 0.01978 0.9482 0.4564
## Residual 47 0.76925 0.98022
## Total 48 0.78477 1.00000
permdisp(ps.epileptic, 'Pheno.Y.N','bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.02029 0.020293 1.5629 10000 0.2151
## Residuals 47 0.61028 0.012985
permdisp(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2646 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0004601 0.00046008 4.9251 10000 0.0292 *
## Residuals 47 0.0043905 0.00009341
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# within epileptic dogs
table(sample_data(ps)$Seizure.Freedom..Y.N., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## No 0 32
## Yes 0 14
table(sample_data(ps)$Seizure.Control..Satisfactory.Unsatisfactory., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## S 0 33
## US 0 13
Does health condition related to dog’s sex? Here we conduct a Chi-squared test.
# chisq test on epi male dog v.s. female.
# not within health condition
tb <- table(sample_data(ps)$Sex, sample_data(ps)$Epileptic.or.Control)
tb
##
## Control Epileptic
## F 33 25
## M 16 24
chisq.test(tb)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tb
## X-squared = 2.0698, df = 1, p-value = 0.1502
plot_ord(ps, 'Sex','MDS','bray')
plot_ord(ps, 'Sex','MDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1364 -- in the phylogenetic tree in the data you provided.
plot_ord(ps, 'Sex','NMDS','bray')
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.242098
## Run 1 stress 0.2453813
## Run 2 stress 0.2402009
## ... New best solution
## ... Procrustes: rmse 0.06907647 max resid 0.2808975
## Run 3 stress 0.2497949
## Run 4 stress 0.2425178
## Run 5 stress 0.2450577
## Run 6 stress 0.2467007
## Run 7 stress 0.2470747
## Run 8 stress 0.2524667
## Run 9 stress 0.2505825
## Run 10 stress 0.2536592
## Run 11 stress 0.2438436
## Run 12 stress 0.2468012
## Run 13 stress 0.2552667
## Run 14 stress 0.2526795
## Run 15 stress 0.2630919
## Run 16 stress 0.2425955
## Run 17 stress 0.2442017
## Run 18 stress 0.2600297
## Run 19 stress 0.2606059
## Run 20 stress 0.244812
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
plot_ord(ps, 'Sex','NMDS','wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2524 -- in the phylogenetic tree in the data you provided.
## Run 0 stress 0.1904303
## Run 1 stress 0.1969604
## Run 2 stress 0.1890081
## ... New best solution
## ... Procrustes: rmse 0.03639876 max resid 0.1910498
## Run 3 stress 0.18442
## ... New best solution
## ... Procrustes: rmse 0.06855117 max resid 0.3274062
## Run 4 stress 0.1843571
## ... New best solution
## ... Procrustes: rmse 0.01999591 max resid 0.1783847
## Run 5 stress 0.1911162
## Run 6 stress 0.1856643
## Run 7 stress 0.1856803
## Run 8 stress 0.1980304
## Run 9 stress 0.1902333
## Run 10 stress 0.1901071
## Run 11 stress 0.1886344
## Run 12 stress 0.19008
## Run 13 stress 0.1863481
## Run 14 stress 0.1856805
## Run 15 stress 0.1860512
## Run 16 stress 0.1904701
## Run 17 stress 0.1898851
## Run 18 stress 0.1925002
## Run 19 stress 0.1863481
## Run 20 stress 0.1907652
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.prop, 'Household + Sex', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.2373 9.999e-05 ***
## Sex 1 0.0589 0.00332 0.5174 0.9774
## Residual 48 5.4604 0.30787
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1628 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 1.21161 0.67613 2.0949 9.999e-05 ***
## Sex 1 0.00201 0.00112 0.1669 0.9978
## Residual 48 0.57836 0.32275
## Total 97 1.79198 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.prop, 'Sex', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01829 0.018289 1.2389 10000 0.2721
## Residuals 96 1.41725 0.014763
permdisp(ps.prop, 'Sex', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV1999 -- in the phylogenetic tree in the data you provided.
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.000131 0.00013143 0.0899 10000 0.7675
## Residuals 96 0.140320 0.00146166
sam$Age..months. <- as.numeric(sam$Age..months.)
ggplot(sam) +
geom_line(aes(x = as.numeric(Household), y = Age..months./12, group = Household)) +
geom_point(aes(x = as.numeric(Household), y = Age..months./12, colour = Epileptic.or.Control)) +
xlab('Household') + ylab('Age in Year')
sam.epi <- sam %>% filter(Epileptic.or.Control == 'Epileptic')
sam.ctl <- sam %>% filter(Epileptic.or.Control == 'Control')
age.diff <- data.frame(Household = sam.epi$Household, Age.Diff = sam.epi$Age..months. - sam.ctl$Age..months.)
ggplot(age.diff) +
geom_bar(aes(x = as.numeric(Household), y = Age.Diff/12), stat = 'identity') +
xlab('Household') + ylab('Age Difference in Year')
permanova(ps.prop, 'Household + Age..months.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2165 0.68881 2.3084 9.999e-05 ***
## Age..months. 19 2.2116 0.12470 1.0557 0.315
## Residual 30 3.3076 0.18650
## Total 97 17.7358 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.prop, 'Household + Age..months.', 'wunifrac')
## Warning in UniFrac(physeq, weighted = TRUE, ...): Randomly assigning root as --
## ASV2330 -- in the phylogenetic tree in the data you provided.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 0.054452 0.67232 2.1861 9.999e-05 ***
## Age..months. 19 0.010971 0.13546 1.1127 0.2822
## Residual 30 0.015567 0.19221
## Total 97 0.080990 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1